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The dynamics of neural networks is often characterized by collective behavior and quasi-synchronous 
events, where a large fraction of neurons fire in short time intervals, separated by uncorrected firing 
activity. These global temporal signals are crucial for brain functioning. They strongly depend on the 
topology of the network and on the fluctuations of the connectivity. We propose a heterogeneous mean- field 
approach to neural dynamics on random networks, that explicitly preserves the disorder in the topology at 
growing network sizes, and leads to a set of self-consistent equations. Within this approach, we provide an 
effective description of microscopic and large scale temporal signals in a leaky integrate-and-fire model with 
short term plasticity, where quasi-synchronous events arise. Our equations provide a clear analytical picture 
of the dynamics, evidencing the contributions of both periodic (locked) and aperiodic (unlocked) neurons to 
the measurable average signal. In particular, we formulate and solve a global inverse problem of 
reconstructing the in-degree distribution from the knowledge of the average activity field. Our method is 
very general and applies to a large class of dynamical models on dense random networks. 



Topology has a strong influence on phases of dynamical models defined on a network. Recently, this topic has 
attracted the interest of both theoreticians and applied scientists in many different fields, ranging from 
physics, to biology and social sciences'"''. Research has focused in two main directions. The direct problem 
aims at predicting the dynamical properties of a network from its topological parameters'. The inverse problem is 
addressed to the reconstruction of the network topological features from dynamic time series" The latter 
approach is particularly interesting when the direct investigation of the network is impossible or very hard to be 
performed. 

Neural networks are typical examples of such a situation. In local approaches to inverse problems'* the 
network is reconstructed through the knowledge of long time series of single neuron dynamics, a methods that 
applies efficiently to small systems only. Actually, the signals emerging during neural time evolution are often 
records of the average synaptic activity from large regions of the cerebral cortex - a kind of observable likely much 
easier to be measured than signals coming from single neuron activities'^ ". Inferring the topological properties of 
the network from global signals is still an open and central problem in neurophysiology. In this paper we 
investigate the possibility of formulating and solving such si global version of the inverse problem, reconstructing 
the network topology that has generated a given global (i.e. average) synaptic-activity field. The solution of such 
an inverse problem could also imply the possibility of engineering a network able to produce a specific average 
signal. 

As an example of neural network dynamics, we focus on a system of leaky integrate-and-fire (LIF) excitatory 
neurons, interacting via a synaptic current regulated by the short-term plasticity mechanism''*''^. This model is 
able to reproduce synchronization patterns observed in in vitro experiments"""*. As a model for the underlying 
topology we consider randomly uncorrelated diluted networks made of N nodes. In general N is considered quite 
a large number, as is the number of connections between pairs of neurons. This suggests that the right framework 
for understanding large-population neural networks should be a mean-field approach, where the thermodyn- 
amic limit, is expected to provide the basic ingredients for an analytic treatment. On the other hand, the 
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way such a thermodynamic Umit is performed may wipe out any 
relation with the topological features that are responsible, for finite N, 
of relevant dynamical properties. 

In Erdos-Renyi directed networks, where each neuron is ran- 
domly and uniformly connected to a finite fraction of the other 
neurons {massive or dense connectivity), the fluctuations of the 
degree determine a rich dynamical behavior, characterized in par- 
ticular by quasi-synchronous events (QSE). This means that a large 
fraction of neurons fire in a short time interval of a few milliseconds 
(ms), separated by uncorrelated firing activity lasting over some tens 
of ms. Such an interesting behavior is lost in the thermodynamic 
limit, as the fluctuations of the connectivity vanish and the "mean- 
field-like" dynamics reduces to a state of fully synchronized neurons 
(e.g., see"). In order to maintain the QSE phenomenology in the large 
N limit, we can rather consider the sequence of random graphs that 

keep the same specific in-degree distribution P(j^ > where ki = ki/N 

is the fraction of incoming neurons connected to neuron i for any 
finite N, similarly to the configuration models™. This way of per- 
forming the thermodynamic limit preserves the dynamical regime of 
QSE and the difference between synchronous and non-synchronous 
neurons according to their specific in-degree k. By introducing expli- 
citly this N ^ 00 limit in the differential equations of the model, we 
obtain a heterogeneous mean-field (HMF) description, similar to the 
one recently introduced in the context of epidemiological spreading 
on networks^''-^''^^. Such mean-field like or HMF equations can be 
studied analytically by introducing the return maps of the firing 
times. In particular, we fmd that a sufficiently narrow distributions 

of P(j^ is necessary to observe the quasi-synchronous dynamical 

phase, which vanishes on the contrary for broad distributions of k. 

More importantly, these HMF equations allow us to design a 
"global" inverse-problem approach, formulated in terms of an integ- 
ral Fredholm equation of the first kind for the unknown p(jcj". 

Starting from the dynamical signal of the average synaptic-activity 
field, the solution of this equation provides with good accuracy the 

p(j^ of the network that produced it. We test this method for very 

different uncorrelated network topologies, where P(j^^ ranges from 

a Gaussian with one or several peaks, to power law distributions, 
showing its effectiveness even for finite size networks. 

The overall procedure applies to a wide class of network dynamics 
of the type 




where the vector Wj represents the state of the site i, F(w,-, 0) is the 
single site dynamics, g is the coupling strength, G(w,) is a suitable 
coupling function and , is the adjacency matrix of the directed 
uncorrelated dense network, whose entries are equal to 1 if neuron 
j fires to neuron i, and 0 otherwise. 

Results 

The LIP model with short term plasticity. Let us introduce LIE 
models, that describe a network of Nneurons interacting via a synap- 
tic current, regulated by short-term-plasticity with equivalent 
synapses"". In this case the dynamical variable of the neuron / is w, 
= (v„ Xj, yi, Zi) where v, is the rescaled membrane potential and x,-, y„ 
and Zj represent the fractions of synaptic transmitters in the 
recovered, active, and inactive state, respectively (%, + yi + Zi = 1). 
Eq. (1) then specializes to: 

v, = fl-v,-h^^e;,7; (2) 



7,. = '- + uXjSi (3) 



The function Sj{t) is the spike-train produced by neuron j, 
Sj{t)= '^m^{^~^)i'^))' where tj{m) is the time when neuron j 
fires its m-th spike. Notice that we assume the spike to be a (5 
function of time. Whenever the potential v,(t) crosses the threshold 
value Vth = 1, it is reset to Vj = 0, and a spike is sent towards its 
efferent neurons. The mechanism of short-term plasticity, that 
mediates the transmission of the field Sj(t), was introduced in"'^ 
to account for the activation of neurotransmitters in neural 
dynamics mediated by synaptic connections. In particular, when 
neuron ; emits a spike, it releases a fraction of neurotransmitters 
uXi(t) (see the second term in the r.h.s. of Eq. (3)), and the fraction 
of active resources y,(f) is increased. Between consecutive spikes of 
neuron /, the use of active resources determines the exponential 
decrease of yi{t), on a time scale Tj^, thus yielding the increment of 
the fraction of inactive resources z,(0 (see the first term on the r.h.s. 
of Eq. (4)). Simultaneously, while z,(t) decreases (see the second term 
on the r.h.s. of Eq. (4)), the fraction of available resources is recovered 
over a time scale Zj.: in fact, from Eq.s (3) and (4) one readily obtains 
Xi{t) =Zi{t) /T^ — uXj{t)Si{t). We assume that all parameters 
appearing in the above equations are independent of the neuron 
indices, and that each neuron is connected to a macroscopic num- 
ber, 0{N), of pre-synaptic neurons: this is the reason why the 
interaction term is divided by the factor N. In all data hereafter 
reported we have used phenomenological values of the rescaled 
parameters: Tin = 0.2, = 133Tin, a = 13, g = 30 and u = 0.5". 

The choice of the value of the external current, a, is quite import- 
ant for selecting the dynamical regime one is interested to reproduce. 
In fact, for a > v^^ = 1, neurons are in a firing regime, that typically 
gives rise to collective oscillations"'^'''^^'^'. These have been observed 
experimentally in mammalian brains, where such a coherent rythmic 
behavior involves different groups of neurons^"". On the other hand, it 
is also well known that in many cases neurons operate in the presence 
of a subthreshold external current'". In this paper, we aim to present a 
method that works irrespectively of the choice of a. For the sake of 
simplicity, we have decided to describe it for a = 1.3, i.e. in a strong 
firing regime. 

Numerical simulations can be performed effectively by transform- 
ing the set of differential equations (2)-(4) into an event-driven 
jjj^pi9.28,29 Erdos-Renyi random graphs, where each link is 
accepted with probability p, so that the average in-degree (fc) = 
pN, the dynamics has been analyzed in detail". Neurons separate 
spontaneously into two different families: the locked and the 
unlocked ones. The locked neurons determine quasi-synchronous 
events (QSE) and exhibit a periodic dynamics. The unlocked ones 
participate in the uncorrelated firing activity and exhibit a sort of 
irregular evolution. Neurons belong to one of the two families 
according to their in-degree fc,. In this topology, the thermodynamic 
limit can be simply worked out. Unfortunately, this misses all the 
interesting features emerging from the model at finite N. Actually, for 
any finite value of N, the in-degree distribution P{k) is centered 
around (k), with a standard deviation ai^-^N^. The effect of disorder 
is quantified by the ratio Oj^lQi), that vanishes for N— » oo. Hence the 
thermodynamic limit reproduces the naive mean-field like dynamics 
of a fully coupled network, with rescaled coupling g pg, that is 
known to eventually converge to a periodic fully synchronous state". 

The LIP model on random graphs with fixed in-degree distri- 
bution. At variance with the network construction discussed in 
previous sections, uncorrelated random graphs can be defined by 
different protocols, that keep track of the in-degree inhomogeneity 
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Figure 1 | Raster plot representation of the dynamics of a network of JV = 
500 neurons with a Gaussian in-degree distribution P(kj 
^(jij = 0.7, (Tj, = 0.077^ . Black dots signal single firing events of neurons 
at time t. Neurons are naturally ordered along the vertical axis according to 
the values of their in-degree. The global field Y ( t) (red line) is superposed 
to the raster plot for comparison (its actual values are multiplied by the 
factor 10', to make it visible on the vertical scale). 



in the thermodynamic limit. In our construction, we fix the 
normalized in-degree probability distribution P(j^^' so that (TkKk) 
is kept constant for any AP". Accordingly, P(i^ is a normalized 
distribution defined in the interval fcE(0,l] (while the number of 
inputs k e (0, N]). In particular, if P(j^ is a truncated Gaussian 

distribution, the dynamics reproduces the scenario observed in''-" for 
an Erdos-Renyi random graph. In fact, also in this case neurons are 
dynamically distinguished into two families, depending on their in- 
degree. Precisely, neurons with k in between two critical values, k^ 
and kc2~(k^, are locked and determine the QSE: they fire with 

almost the same period, but exhibit different (fc-dependent) phases. 
All the other neurons are unlocked and fire in between QSE 
displaying an aperiodic behavior. Notice that the range 0 < fc < A:^ 
corresponds to the left tail of the truncated Gaussian distribution; 
accordingly, the large majority of unlocked neurons is found in the 
range kc2 <k<l (see Fig. 1). In order to characterize the dynamics at 
increasing N, we consider for each neuron its inter-spike-interval 
(ISI), i.e. the lapse of time in between consecutive firing events. In 

Fig. 2 we show the time-average of ISI vs k, or ISI (jcj . One can clearly 
observe the plateau of locked neurons and the crossover to unlocked 
neurons at the critical values k^ and fcc2- 

Remarkably, networks of different sizes {N = 500, 5000 and 

20000) feature the same ISI (jc^ for locked neurons, and almost the 

same values of k^ and kc2 ■ There is not a sharp transition from locked 
to unlocked neurons, because for finite N the behavior of each neu- 
ron depends not only on its k, but also on neighbor neurons sending 
their inputs. Nevertheless, in the inset, the crossover appears to be 
sharper and sharper for increasing N, as expected for true critical 

points. Furthermore, the fluctuations of ISl(JcJ over different reali- 
zations, by PI ) , of three networks of different size exhibit a peak 



around k^ and kc2, while they decrease with N as ~ N~"^ (data not 
shown). Thus, the qualitative and quantitative features of the QSE at 
finite sizes are expected to persist in the thermodynamic limit, where 



fluctuations vanish and the dynamics of each neuron depends only 
on its in-degree. 

Heterogeneous mean field equations. We can now construct the 
Heterogeneous Mean-Field (HMF) equations for our model by 
combining this thermodynamic limit procedure with a consistency 
relation in Eqs. (2)-(4). The input field received by neuron i is 
Y, = l/Nj2j<^ijyj{t) = '^/NY,,ei{i)yp where I{i) is the set of fc,- 
neurons transmitting to neuron For very large values of N the 
average field generated by presynaptic neurons of neuron i.e. 

1 y ^i^^ '"'"^ approximated with I/N^^.jK;, where the 

second sum runs over all neurons in the network {mean-field 
hypothesis). Accordingly, we have Y,- = (fc,/N)(l//c,) ^^.^ 

ki(^ /N^^^/yj^: as a consequence in the thermodynamic limit the 
dynamics of each neuron depends only on its specific in-degree. In 
this limit, fc, becomes a continuous variable independent of the label 
taking values in the interval (0,1]. Therefore, aU neurons with the 
same in-degree k follow the same dynamical equations and we can 
write the dynamical equations for the class of neurons with in-degree 
k. Finally, we can replace Y,- with kY[t), where 



dk 



The HMF equations, then, read 

h^t)='^-nit)+gkY{t) 



+uii-rk{t)-zi{t))Sj^{t) 



(5) 



(6) 



(7) 



(8) 



where v^y^. and z^. are the membrane potential, fraction of active and 
inactive resources of the class of neurons with in-degree k, 
respectively. Despite this set of equations cannot be solved 
explicitly, they provide a great numerical advantage with respect to 
direct simulations of large systems. Actually, the basic features of the 
dynamics of such systems can be effectively reproduced (modulo 

finite-size corrections) by exploiting a suitable sampling of P(k\ 



For instance, one can sample the continuous variable /ce[0,1] into M 



discrete values kj in such a way that 



P[k\dk is kept fixed 



(importance sampling). Simulations of Equations (5)-(8) show 
that the global field field Y {t) is periodic and the neurons split into 
locked and unlocked. Locked neurons feature a firing time delay with 
respect the peak of Y (t), and this phase shift depends on the in- 
degree k. As to unlocked neurons, that typically have an in-degree 
k>{k\ they follow a complex dynamics with irregular firing times. 



In Fig. 2 (black dots) we compare ISI ykj , obtained from the HMF 

equations for M = 307, with the same quantity computed by direct 
simulations of networks up to size N = 2 X 10''. The agreement is 
remarkable evidencing the numerical effectiveness of the method. 

The direct problem: stability analysis and the synchronization 
transition. In the HMF equations, once the global field Y (f) is 
known, the dynamics of each class of neurons with in-degree k can 
be determined by a straightforward integration, and we can perform 
the stability analysis that Tsodyks et al. applied to a similar modeP". 
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Figure 2 | Time average of inter-spike intervals ISI [kj \s k from a 



Gaussian distribution with (icj =0.7 and ct^ = 0.07/" and for three 
networks with N = 500 (blue triangles), N = 5000 (red diamonds), N = 
20000 (green squares). For each size, the average is taken over 8 different 
realizations of the random network. We have also performed a suitable 
binning over the values of k, thus yielding the numerical estimates of the 
critical values fcci~0.49 and k^2~0-70. In the inset we show a zoom of the 
crossover region close to ^fc^ = 0.7. Black dots are the result of simulations 
of the mean field dynamics (see Eq.s (6)-(8)) with M = 307. 

As an example, we have considered the system studied in Fig. 2. The 
global field Y (t) of the HMF dynamics has been obtained using the 

importance sampling for the distribution P(j^- have fitted Y (f) 

with an analytic function of time Y^t), continuous and periodic in 
time, with period T. Accordingly, Eq. (6) can be approximated by 

viit)='^-^l(t)+gkYf{t). (9) 

Notice that, by construction, the field Yj(t) features peaks at times t = 
nT, where n is an integer. In this way we can represent Eq. (9) as a 
discrete single neuron map. In practice, we integrate Eq.(9) and 
determine the sequence of the (absolute value of the) firing time- 
delay, ( «) , of neurons with in-degree k with respect to the reference 
time nT. The return map Rj. of this quantity reads 

(10) 



In Fig. 3 we plot the return map of the rescaled firing time-delay 
ti{n)/T for different values of k. We observe that in-degrees k 
corresponding to locked neurons (e.g., the brown curve) have two 
fixed points and t^, i.e. two points of intersection of the curve with 
the diagonal. The first one is stable (the derivative of the map Rj^ is 
<1) and the second unstable (the derivative of the map Rj. is >1). 
Clearly, the dynamics converges to the stable fixed point displaying a 
periodic behavior. In particular, the firing times of neurons k are 
phase shifted of a quantity t| with respect the peaks of the fitted 
global field. The orange and violet curves correspond to the dyna- 
mics at the critical in-degrees k^ and kc2 where the fixed points 
disappear (see Fig. (2)). The presence of such fixed points 
influences also the behavior of the unlocked component (e.g., the 
red and light blue curves). In particular, the nearer k is to kd or to kc2, 
the closer is the return map to the bisector of the square, giving rise to 
a dynamics spending longer and longer times in an almost periodic 
firing. Afterwards, unlocked neurons depart from this almost 
periodic regime, thus following an aperiodic behavior. As a 
byproduct, this dynamical analysis allows to estimate the values of 
the critical in-degrees. For the system of Fig. 1, A:ci=0.48 and 




0.4 0.6 
t(n) 



Figure 3 | The return map R^. in Eq. ( 10) of the rescaled variables t^(«)/T 
for different values of fc, corresponding to lines of different colors, 
according to the legend in the inset: the black Kne is the bisector of the 
square. 



fcc2 = 0.698, in very good agreement with the numerical 
simulations (see Fig. 2). 

StiU in the perspective of the direct problem, the HMF equations 
provide further insight on how the network topology influences the 
dynamical behavior. We have found that, in general, the fraction of 

locked neurons increases as P(j^ becomes sharper and sharper, 

while synchronization is eventually lost for broader distributions. 

In Fig. 4 we report the fraction of locked neurons,/; = ^( ^ ) 
as a function of the standard deviation deviation Cj., for different 
kinds of (single or double-peaked Gaussian, power law) in 

the HMF equations. For the truncated power law distribution, we 
set P^fc^ ~ o{k — k^^ fc^". In all cases, there is a critical value of a-^ 
above which // vanishes, i.e. QSE disappear. The generality of this 
scenario points out the importance of the relation between P^fc^ and 
the average synaptic field Y(t). 

The inverse problem. The HMF approach allows to implement the 
inverse problem and leads to the reconstruction of the distribution 

from the knowledge of Y(i). If the global synaptic activity field 

Y(j) is known, each class of neurons of in-degree k evolves according 
to the equations: 



yi{i)=a-y^^{t)+gkY{t) 



+ M(l-y^(0 -2^(0)5^(0 



n(0 



(11) 



(12) 



(13) 



Notice that the variable v{t), j(t), z{t) can take values that differ from 
the variables generating the field Y{t), i.e. v{t), y{t), z{t), as they start 
from different initial conditions. However, the self consistent relation 
for the global field Y{t) implies: 
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Figure 4 | The fraction of locked neurons, fi, as a function of the standard 
deviation (7^ of the distributions: truncated Gaussian with (jcj =0.7 
(black dots); truncated superposition of two Gaussians (both with 
standard deviation 0.03), one centered at fci = 0.5 and the other one at a 
varying value ^2, that determines the overall standard deviation (blue 
squares); truncated power law distribution with fcmi„ = 0.1 (red 
diamonds). In the last case the value of the standard deviation is changed 
by varying the exponent a, while the average (jcj changes accordingly. 
Lines have been drawn to guide the eyes. 



P{k)Y,it)dk. 



(14) 



If Y{t) and Yi{t) are known, this is a Fredholm equation of the first 

kind in P^A:^^'. In the general case of Eq. (1), calling £(t) the global 

measured external field, the evolution equations corresponding to 
Eq.s (11)-(13) read 



and the Fredholm equation for the inverse problem is 



E(t) 



P k 



:)G(w^(f))dfc. 



15 



(16) 



In the case of our LIF model, as soon as a locked component exists, 
Eq. (14) can be solved by a functional Montecarlo minimization 

procedure applied to a sampled P(j^- variance with the direct 
problem, P (Jc^ is the unknown function and, accordingly, we have to 
adopt a uniform sampling of the support of k. A sufficiently fine 
sampling has to be used for a confident reconstruction of p(^k 
(See Methods section). 

To check our inverse method, we choose a distribution P{^k 

evolve the system and extract the global synaptic field Y(t). We then 
verify if the procedure reconstructs correctly the original distribution 

p(j^ ■ In panels (a), (b) and (c) of Fig. 5 we show examples in which 
Y(f) has been obtained from the simulation of the HMF with differ- 
ent P^fc^ (Gaussian, double peak Gaussian and power law). We can 
see that the method determines confidently the original distribution 
p(j^ ■ Notice that the method fails as soon as the locked component 

disappears, as explained in the methods section. Remarkably, the 
method can recognize the discontinuity of the distribution in 



P(K)8l 



P(k) 




0.2 0, 



0.2 0.4 0.6 0.8 
1 



Figure 5 | Inverse problem for Pl^kj from the global field Yit). 

Panels (a), (b) and (c) show three distributions of the kind considered in 
Fig. (4) (black continuous curves) for the HMF equations and their 
reconstructions (circles) by the inverse method. The parameters of the 
three distributions are = 0.043, k2 = 0.7 and a = 4.9. In panel (d) we 



show the reconstruction (crosses) o{ Pykj (black continuous line) by the 
average field Y{ t) generated by the dynamics o{n finite size network with N 
= 500. 

k = k^in and the value of the exponent of the power law a = 4.9. 
Finally, in panel (d) of Fig. 5, we show the result of the inverse 

problem for the distribution P(k^ obtained from a global signal 

generated by a finite size realization with N = 500 and (k) = 350. 
The significant agreement indicates that the HMF and its inverse 
problem are able to infer the in-degree probability distribution 

P (ji^ even for a realistic finite size network This last result is par- 
ticularly important, as it opens new perspectives for experimental 
data analysis, where the average neural activity is typically measured 
from finite size samples with finite but large connectivity. 

Discussion 

The direct and inverse problem for neural dynamics on random 
networks are both accessible through the HMF approach proposed 
in this paper. The mean-field equations provide a semi-analytic form 
for the return map of the firing times of neurons and they allow to 
evidence the effects of the in-degree distribution on the synchron- 
ization transition. This phenomenology is not limited to the LIF 
model analyzed here and it is observed in several neural models on 
random networks. We expect that the HMF equations could shed 
light on the nature of this interesting, but still not well understood, 
class of synchronization transitions^'''^^'^'. The mean field nature of 
the approach does not represent a limitation, since the equations are 
expected to give a faithful description of the dynamics also in net- 
works wdth finite but large average in-degree, corresponding to the 
experimental situation observed in many cortical regions''''. 

The inverse approach, although tested here only on numerical 
experiments, gives excellent results on the reconstruction of a wide 
class of in-degree distributions and it could open new perspectives on 
data analysis, allowing to reconstruct the main topological features of 
the neural network producing the QSE. The inversion appears to be 
stable with respect to noise, as clearly shown by the effectiveness of 
the procedure tested on a finite realization, where the temporal sig- 
nals of the average synaptic activity is noisy. We also expect our 
inverse approach to be robust with respect to the addition of limited 
randomness in leaky currents, and also with respect to a noise com- 
patible with the signal detection from instruments. Finally, the 
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method is very general and it can be applied to a wide class of 
dynamical models on networks, as those described in Eq. (1). 

Further developments will allow to get closer to real experimental 
situations. We mention the most important, i.e. the introduction of 
inhibitory neurons and the extension of our approach by taking into 
account networks with degree correlations', that are known to char- 
acterize real structures, and sparse networks. The HMF equations 
appear to be sufficiently flexible and simple to allow for these 
extensions. 

Methods 

Montecarlo approach to the inverse problem. In this section we provide details of 
the algorithmic procedure adopted for solving the inverse problem, i.e. reconstructing 

the distribution p(jcj from Eq. (14). In the HMF formulation, the field Y(t) is 

generated by an infinite number of neurons and /: is a continuous variable in the 
interval (0, 1]. In practice, we can sample uniformly this unit interval by L disjoint 
subintervals of length 1/L, labelled by the integer i. This corresponds to an effective 
neural index i, that identifies the class of neurons with in-degree ki — i/L. In this way 
we obtain a discretized definition converging to Eq.(14) for L ^\ 

nt)= {p(fc)yiWd^- jEJ'(^*)y;i,(t)- (17) 

In order to improve the stability and the convergence of the algorithm by smoothing 
the fluctuations of the fields y^ (0' convenient to consider a coarse-graining of 
the sampling by approximating 7(f) as follows 



(18) 



where ^y^,(Oy average oiLIV synaptic fields of connectivity fcG|^fc,-,fc, _|_i J . This 
is the discretized Fredholm equation that one can solve to obtain ^{^i^ from the 
knowledge of {yi\t)^ and Y{t). For this aim we use a Monte Carlo (MC) 
minimization procedure, by introducing at each MC step, n, a trial solution, P„ {li-^ , 
in the form of a normalized non-negative in-degree distribution. Then, we evaluate 



the field Y„{i) and the distance 7„ defined as: 

L'-l 

L' 



Y„ t,P„ k, 



:))-y{t))] 



-dt. 



(19) 



(20) 



The time interval [f i, has to be taken large enough to obtain a reliable estimate of 
For instance, in the case shown in Fig. 1, where Y{t) exhibits an almost periodic 
evolution of period T ~ 1 in the adimensional units of the model, we have used t2 — f i 
= 10. The overall configuration of the synaptic fields, at iteration step n + 1, is 
obtained by choosing randomly two values kj and fcj, and by defining a new trial 
solution Pn + \{k^ ^Pn ~^^^\k ~^^\~kr provided both + i ^fc^j and 

Pn + i (j^i^ are non-negative, we increase and decrease P„ (jcj^ of the same amount, e, 
in kj and kj respectively. A suitable choice is O(l0~*). Then, we evaluate 
7«+i(p„+i(fc,)): If7„ + i (^,, + 1 (^f)) <7„(p„(fci)) the step is accepted i.e. 
P„_i_i — + otherwise P^+i — Pn- This MC procedure amounts to the imple- 
mentation of a zero temperature dynamics, where the cost function 7„ {^n ^ can 
only decrease. In principle, the inverse problem in the form of Eq.(18) is solved, i.e. 

Yn [t^Pn (^/)) ^Y{t), if7„ {Pn (^;) ) 

— 0. In practice, the approximations introduced 
by the coarse- graining procedure do not allow for a fast convergence to the exact 
solution, but Pf, {ki^ can be considered a reliable reconstruction of the actual 
already for 7« < 10"^. We have checked that the results of the MC procedure are quite 
stable with respect to different choices of the initial conditions Pq {ki^ , thus con- 
firming the robustness of the method. We give in conclusion some comments on the 
very definition of the coarse-grained synaptic field ^y^, ( f ) ^ Since small differences in 
the values of ki reflect in small differences in the dynamics, for not too large intervals 
|^fc,-,/c,_]_ 1 j the quantity ^y^ ( t) ^ can be considered as an average over different initial 
conditions. For locked neurons the convergence of the average procedure defining 



\yi.{t)j is quite fast, since all the initial conditions tend to the stable fixed point, 

identified by the return map described in the previous subsection. On the other hand, 
the convergence of the same quantity for unlocked neurons should require an average 
over a huge number of initial conditions. For this reason, the broader is the distri- 
bution, i.e. the bigger is the unlocked component (see Fig. 4), the more computa- 
tionally expensive is the solution of the inverse problem. This numerical drawback for 
broad distributions emerges in our tests of the inversion procedure described in Fig. 5. 
Moreover, such tests show that the procedure works insofar the QSE are not neg- 
ligible, but it fails in the absence of the locking mechanism. In this case, indeed, the 

global field Y{t) is constant and also ^y^ [t) ^ become constant, when averaging over a 

sufficiently large number of samples. This situation makes Eq.(18) trivial and useless 

to evaluate ^{^i^ ■ We want to observe that, while in general y^ (f) ^^y^. (t), one can 

reasonably expect that ^y^ (f)^ ^ ^^H^ good approximation of (yiX^)"^- This 
remark points out the conceptual importance of the HMF formulation for the 
possibility of solving the inverse problem. 
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